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Abstract 

Computer  modeling  is  a  key  tool  in  understanding  wave  phenom¬ 
ena  in  nature.  It  is  used  for  instance  in  seismology  for  exploring  the 
consequences  of  the  Earth’s  fine  structure  and  in  electromagnetism 
for  modeling  of  radar  and  microwave  signals.  Iterative,  model  driven 
estimation  of  material  structure  and  constitutive  parameters  is  an  oft- 
proposed  application  of  computer  modeling,  but  requires  very  fast 
simulations  to  achieve  its  potential.  These  simulations  are  computa¬ 
tionally  challenging  especially  in  three  dimensions. 

This  project  studied  a  variety  of  approaches  to  improving  the  ef¬ 
ficiency  and  fidelity  of  wave  simulations.  It  also  developed  a  basic 
software  approach  to  material  parameter  estimation.  This  software 
is  designed,  for  example,  to  assist  in  the  interpretation  of  Ground 
Penetrating  Radar  data  for  civil  and  environmental  engineering  appli¬ 
cations. 
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Research  Topics 

Computer  modeling  is  a  key  tool  in  understanding  wave  phenomena  in  na¬ 
ture.  It  is  used  for  instance  in  seismology  for  exploring  the  consequences 
of  the  Earth’s  fine  structure  and  in  electromagnetism  for  modeling  of  radar 
and  microwave  signals.  Iterative,  model  driven  estimation  of  material  struc¬ 
ture  and  constitutive  parameters  is  an  oft-proposed  application  of  computer 
modeling,  but  requires  very  fast  simulations  to  achieve  its  potential.  These 
simulations  are  computationally  challenging  especially  in  three  dimensions. 

This  project  began  with  the  study  of  a  particular  simulation  approach, 
which  we  called  the  Wavelet/Galerkin  method.  The  use  of  wavelets  as  ba¬ 
sis  function  in  a  Galerkin  method  leads  to  explicit  time  stepping  schemes, 
combining  the  computational  efficency  of  finite  difference  methods  with  the 
analytical  simplicity  of  finite  element  methods.  We  aimed  to  gain  insight 
into  the  behaviour  of  these  methods,  especially  in  comparison  with  standard 
finite  element  and  finite  difference  methods  for  the  same  problems. 

The  application  of  wave  modeling  to  material  property  estimation  formed 
the  other  focus  of  this  project,  especially  the  estimation  of  electromagnetic 
properties  of  soil  and  rock  for  civil  and  environmental  engineering  projects, 
using  Ground  Penetrating  Radar  (GPR).  The  proposed  approach  differed 
from  other  reported  efforts  in  this  topic  in  the  use  of  a  full-blown  numerical 
simulator  for  radar  waves,  driving  optimization  software  for  data  fitting. 


Progress  in  Numerical  Wave  Simulation 

Work  complete  before  the  beginning  of  the  project  determined  that  at  least 
a  considerable  subset  of  wavelet-Galerkin  methods  are  actually  equivalent  to 
explicit  finite  difference  methods,  though  in  some  cases  with  novel  choices  of 
stencils,  especially  for  order  higher  than  four.  For  the  application  to  GPR, 
we  confined  our  attention  for  the  duration  of  this  project  to  schemes  of  order 
four,  hence  to  standard  finite  difference  methods.  Several  interesting  ques¬ 
tions  about  wavelet-Galerkin  methods  remain  open,  and  were  not  addressed 
in  our  work,  notably  whether  these  methods  result  in  a  superior  treatment 
of  boundary  conditions  for  systems,  such  as  free  surfaces  in  elastodynamics. 

Computation  of  wave  phase  without  other  details  of  wave  motion  is  pos¬ 
sible  through  solution  of  the  eikonal  equation  from  geometric  optics.  Finite 
difference  methods  developed  in  the  Pi’s  group  and  elsewhere  perform  these 
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calculations  with  unprecedented  efficiency.  S.  Kim  enhanced  this  efficiency 
considerably  by  observing  a  superconvergence  property.  Some  approaches  to 
material  property  estimation  requires  the  simulation  of  the  perturbational 
eikonal  equation  and  its  adjoint  problem.  In  carrying  out  these  perturba¬ 
tional  and  adjoint  calculations,  A.  Sei  discovered  a  disparity  between  numer¬ 
ical  accuracy  of  the  perturbational  simulator  and  the  adjoint  perturbational 
simulator,  which  turns  out  to  be  fundamental  and  unavoidable. 

While  this  project  emphasizes  broadband  GPR,  frequency  domain  sim¬ 
ulation  and  analysis  remain  relevant  for  various  reasons.  S.  Kim  developed 
further  work  begun  in  his  thesis  (under  Prof.  Jim  Douglas  at  Purdue)  on 
domain  decomposition  for  the  Helmholtz  equation.  Finite  element  and  finite 
difference  discretizations  of  this  equation,  a  prototype  for  frequency-domain 
Maxwell,  yield  very  large  nonsymmetric  illconditioned  linear  systems,  for 
which  standard  iterative  methods  converge  slowly  or  fail  to  converge  alto¬ 
gether.  B.  Despres  pioneered  the  mathematics  of  domain  decomposition  to 
produce  preconditioners  for  these  systems.  Kim  has  answered  many  questions 
about  the  discrete  aspects  of  domain  decomposition  left  open  by  Despres. 
He  has  produced  a  class  of  methods  with  optimal  convergence  properties, 
i.e.  depending  weakly  or  not  at  all  on  discretization,  by  combining  domain 
decomposition,  multigrid,  and  various  Krylov  iterations. 

Interfaces  in  Earth  materials  exist,  across  which  material  parameters  take 
large  jumps.  An  example  is  the  water  table,  at  which  the  speed  of  compres- 
sional  seismic  waves  increases  by  a  factor  of  up  to  five,  and  the  electrical 
conductivity  jumps  by  orders  of  magnitude  in  some  cases.  Finite  difference 
methods  are  generally  subject  to  degradation  of  accuracy  in  the  presence  of 
parameter  discontinuities  such  as  these.  C.  Zhang,  continuing  his  thesis  work, 
devised  several  increasingly  sophisticated  variants  of  the  immersed  interface 
method ,  which  maintains  the  accuracy  of  finite  difference  methods  even  in 
the  presence  of  jump  discontinuities  in  parameters.  This  approach  is  being 
further  developed  by  others,  and  will  be  very  important  in  realistic  modeling 
of  Earth  response  and  estimation  of  parameters  from  field  measurements. 

To  initialize  the  work  on  GPR  imaging,  Zhang  also  produced  an  overview 
report  summarizing  the  state  of  the  art  in  GPR  simulation. 

The  work  overviewed  in  this  section  is  reported  in  [13,  15,  14,  3,  8,  6,  7, 
9,  20,  18,  17,  16,  21,  19] 
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Progress  in  Earth  parameter  estimation 

Perturbational  and  adjoint  perturbational  computations  for  wave  equations 
are  the  drivers  of  model-based  parameter  estimation,  as  envisioned  for  GPR 
in  the  second  goal  of  this  project.  The  construction  of  these  simulators  is 
complicated:  several  related  code  modules  must  function  in  consort.  Espe¬ 
cially  production  of  the  adjoint  code  is  tedious  and  error  prone.  A  number 
of  researchers  in  various  areas  (for  example  oceanography  and  meteorology, 
where  many  such  simulator-driven  parameter  estimation  problems  arise)  have 
developed  so  called  Automatic  Differentiation  (“AD”)  software,  which  auto¬ 
mates  the  production  of  perturbational  and  adjoint  perturbational  code  from 
basic  simulator  code.  Certain  technical  problems  render  this  automation  ex¬ 
ceedingly  difficult  for  whole  simulators,  whereas  most  of  the  packages  can 
handle  the  submodules  which  implement  single  steps.  This  project  saw  the 
development  of  two  generations  of  template  code,  which  fully  implements 
an  efficient  strategy  for  perturbational  and  adjoint  code  for  time  stepping 
methods,  except  for  the  single  step  code,  the  latter  being  supplied  through 
use  of  AD  tools.  J.  Blanch  constructed  the  first,  all-Fortran  template,  and 
demonstrated  its  use  in  the  context  of  viscoacoustic  simulation  and  parameter 
estimation.  C.  Zhang  and  the  PI  collaborated  to  develop  a  second  generation 
object  oriented  template  in  C++,  which  we  successfully  coupled  to  Fortran 
single  step  code  produced  by  AD  tools. 

The  second  generation  template  actually  produces  class  definitions  of  op¬ 
erator  objects  as  defined  by  the  Hilbert  Class  Library  (“HCL”),  a  C++ 
library  for  linking  complex  simulation  with  optimization  code.  This  library 
(development  of  which  is  supported  by  continuing  grants  from  NSF  and  DoE) 
is  designed  to  overcome  a  difficulty  often  faced  by  scientific  programmers,  in 
particular  by  the  investigators  in  this  project.  Procedural  library  code  for 
optimization,  linear  algebra,  etc.,  however  high  in  quality,  is  difficult  to  link 
to  complex  simulators,  because  the  latter  involve  data  structures  incompati¬ 
ble  with  the  primitive  types  used  by  the  library  code.  HCL  overcomes  these 
difficulties  by  separating  cleanly  the  levels  of  abstraction  required  by  the  two 
types  of  task  (simulation  vs.  optimization),  using  object  oriented  program¬ 
ming  techniques. 

Zhang  completed  this  construction  for  TE  polarized  Maxwell’s  equations. 
Unfortunately  the  project  ended  before  extensive  testing  of  the  parameter  es¬ 
timation  capabilities  of  the  package  could  take  place.  Another  very  important 
part  of  the  original  agenda  which  did  not  receive  any  useful  attention  in  the 
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project  was  antenna  modeling.  Both  of  these  tasks  would  need  prosecution 
before  any  serious  evaluation  of  the  proposed  approach  can  be  carried  out. 
The  work  overviewed  in  this  section  is  reported  in  [1,  2,  10, 12,  3,  11,  5,  4] 
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